Geostatistics is a discipline concerned with spatially distributed random variables (also called “regionalized variables”), usually applied to problems in the earth sciences, such as estimation of mineral reserves and delineation of mineral deposits, hydrocarbon reservoirs, and groundwater aquifers. Typically it makes use of two-point statistics summarized in a variogram. Multipoint (or multiple-point) geostatistics (MPS) differs from the rest of variogram-based geostatistics primarily in that it characterizes spatial variability using patterns (sets of points and their configurations) that involve higher order (much greater than order 2) statistics.
One of the goals of multipoint geostatistics is simulation, namely the generation of numerical values along a line, on a surface, or in a volume, such that the set of values match certain given spatial correlation or structural properties (usually derived from a data set called an “analog” or “training image” or “template”) while optionally (in the case called “conditional simulation”) matching predetermined data. In practice, the “analog” may be, for example, a well-known rock volume that is statistically similar to a currently uncharacterized oil reservoir being delineated, and the predetermined data to be matched may be lithology observations at wells, or probabilities of lithologies derived from seismic data. MPS simulations are developed to match two types of constraints: absolute constraints are matched exactly; partial constraints are matched probabilistically, as closely as possible to the constraint value, but they are not guaranteed to match exactly. In the case of petroleum reservoir modeling, examples of absolute constraints are typically data acquired in wells or geological outcrops. Partial constraints are typically derived from seismic data; 1D, 2D, or 3D interpreted spatial trend information; geological category probability fields; or rotation and affinity (or scale) constraints on the size and direction of geological features. Such data are used in a stochastic modeling process to generate one-dimensional (1D), two-dimensional (2D), and/or three-dimensional (3D) spatial distribution of geological categories or rock properties. Since there is a random component involved in MPS simulations, individual stochastic realizations of property fields created by MPS algorithms differ, but the ensemble of realizations provide modelers with improved quantitative estimates of the spatial distribution and uncertainty of values in a modeled volume of interest.
Multipoint geostatistical methods have been demonstrated to be computationally feasible and have been tested on real datasets as set forth in i) Strebelle, “Conditional simulation of complex geological structures using multiple-point statistics”, Mathematical Geology, v. 34, n. 1, 2002, pp. 1-22, ii) Strebelle et al., “Modeling of a deepwater turbidite reservoir conditional to seismic data using principal component analysis and multiple-point geostatistics,” SPE Journal, Vol. 8, No. 3, 2003, pp. 227-235, and iii) Liu et al., “Multiple-point simulation integrating wells, three-dimensional seismic data, and geology,” American Association of Petroleum Geologists Bulletin v. 88, no. 7, 2004, pp. 905-921.
Multipoint geostatistical methods use a numerical training image to represent the spatial variability of geological information. The training image provides a conceptual quantitative description of the subsurface geological heterogeneity, containing possibly complex multipoint patterns of geological heterogeneity. Multipoint statistics conditional simulation anchors these patterns to well data (and/or outcrop data) and to the seismic-derived information (and/or probability field information or constraint grid(s)). An example of such method is described in US-2007-0014435, assigned to Schlumberger Technology Corporation.
Geostatistics relies on the well-known concept of random variables. In simple terms, continuous or discrete properties at various spatial locations are largely unknown or uncertain; hence each property of interest at each spatial location is figured into a random variable whose variability is described by a probability function. In order to perform any type of geostatistical simulation, one requires a decision or assumption of stationarity. In multipoint geostatistical methods, the use of training images is bound by the principle of stationarity as described by Caers, J., and T. Zhang, 2004, “Multiple-point geostatistics: a quantitative vehicle for integrating geologic analogs into multiple reservoir models”, in M. Grammer, P. M. Harris, and G. P. Eberli, eds., Integration of Outcrop and Modern Analogs in Reservoir Modeling, Memoir 80: Tulsa, Okla., AAPG. A random spatial field is said to be stationary if all of its statistical parameters are independent of location (invariant according to any translation). In the case of training images, this stationarity can consist of, but is not limited to, orientation stationarity, where directional elements do not rotate across the training image; and scale stationarity (where the size of elements on the image does not change across the training image).
One multipoint geostatistics method is well known in academia and industry by the name of “Single Normal Equation Simulation” (SNESIM) (Strebelle, S., 2000, “Sequential simulation drawing structures from training images, PhD thesis, Stanford University, 200p). The SNESIM method is generally considered useful for practical applications such as modeling categorical or discrete data types, especially for categorical data in 3D property modeling. In the SNESIM method, the conditional probability density function of all categories at one point is computed using knowledge of the value at a number of nearby points and statistics provided by the training image. SNESIM works with discrete values only (i.e., a finite and usually small number of categories, such as for example five different rock types). Assume there are two categories to be simulated: A (“non-channel”) and B (“channel”). The training image contains a complete representation (i.e., an example) of a spatial distribution of A and B (for example, see FIG. 2). Assume further that the category present at 4 points (points 1, 2, 3, and 4 in FIG. 3), which could represent wells, is known. Conceptually, the SNESIM method computes the probability of finding categories A or B at an unsampled point (point u in FIG. 3) by scanning the training image for all occurrences of the “pattern” (points 1, 2, 3, and 4 in FIG. 3) (that is, the pattern is defined by the relative spatial locations of all five points, and by the categories found at the four points where a value already has been assigned). If five such pattern occurrences (called replicates) are found, and 4 out of 5 replicates show category B at the relative position of the unknown point (u), then the method computes the probability of B at the unknown point, given the particular five-point pattern, to be 4/5 or 80% (while that of A is set to 1/5 or 20%). Furthermore, the method can assign a category to the unknown point by randomly drawing a value from a distribution with 80% probability of B and 20% probability of A if there is no any replicates of the four-point pattern or its sub-patterns found in the training image.
In practice, the SNESIM uses pixel-based (or voxel-based in 3D) sequential simulation method. It starts with a volume to be modeled where all property values are unassigned at all grid cell locations, or one that contains only a few data points to be matched. These volumes are usually represented by a Cartesian grid, where each grid location is called a cell. First, SNESIM decides on a random path for visiting each unassigned cell once and only once. In the first cell, the method searches for nearby known points within a search volume, usually an ellipsoid or rectangular volume around the known point. If one or more known (or already assigned) cells are found, it proceeds as to the next point on the random path as described above to find the probability of finding categories A or B at the unknown point. Armed with the probabilities for each category at the unknown point, the method randomly draws a value (category A or B) from the inferred probability distribution and assigns it to the unknown point. The process repeats at the next cell in the initially-assigned random path until all cells are assigned.
Such methodology was well known in the early 1990's (before it was known as “SNESIM”) (Guardiano, F., and R. M. Srivastava, 1993, Multivariate geostatistics: beyond bivariate moments, in A. Soares, ed., Geostatistics-Troia, v. 1: Dordrecht, Netherlands, Kluwer Academic Publications, p. 133-144). One of the limitations of the first MPS approach, however, was that it was extremely computationally intensive to consult the training image multiple times. Precisely speaking, the training image has to be scanned per simulation node; hence the original SNESIM algorithm proposed by Guardiano and Srivastava was very CPU intensive. In 2000, Strebelle developed a technique to store the information contained in the training image in a special tree-like structure that reduced computations enormously (Strebelle, S., 2000, Sequential simulation drawing structure from training images.: PhD Thesis, Stanford University, Stanford, Calif., USA). With this improvement, the methodology was commonly referred to as the SNESIM method, an algorithm that actually traded CPU with memory. Therefore, coupled with the steadily increasing power of digital computers, SNESIM brought the image-consultation problem down to reasonable but still significant processing times and computer RAM memory requirements.
The SNESIM method is based on a tree structure or “search tree” that stores, in a computer's memory, the spatial geometric relationships between values of a discrete variable to be modeled. In the SNESIM method, the search tree is handled as a typical generic tree data structure containing a number of nodes each pointing to one or more other nodes in a tree-like arrangement (FIG. 4). However, building and consulting this structure is a bottleneck of software efficiency, because this structure tends to be very large to contain enough information for modeling property distributions in several dimensions adequately. This is especially the case if the search tree is built for a training image containing more than 2 categorical variables in 3-dimensional space, as is common in real-world modeling applications. To achieve better reproduction of complex patterns from a training image, a large local window is required by SNESIM, to scan the training image. This typically leads to a big search tree in 3D and for training images with multiple facies. Although the method allows the user to set parameters that control and improve this continuity, setting them to the values that would be needed for adequate three (or more)-dimensional modeling, causes the method either to take an impractical amount of time to run, to fail due to lack of RAM memory even on large modern computers, or to produce non-geological modeling artifacts (particularly when modeling in three or more dimensions). Consequently, the resulting model simulation outputs (realizations) produced using the currently known implementation of SNESIM commonly exhibits unacceptable continuity of spatial features.
Furthermore, the SNESIM method shows a significant computational drawback, because it has not been realized that in practice the SNESIM search tree structure does not really resemble a tree, but rather consists of many long non-branching sequences (or “runlengths”) (see FIG. 6, left side). Handling these in the same manner as a conventional tree structure is very inefficient.
Additionally, the SNESIM method exhibits an additional computational drawback when attempting to simulate non-stationary spatial trends such as rotational and scaling transformations. To mimic local affine transformation information (such as the effects of changes in geological sedimentation directions or the spatial changes in geological feature aspect ratios), SNESIM as published needs to divide the simulation grid into different regions of constant scaling and/or rotational transformation. The method then creates a new search tree for each affine transformation region. For common real-world applications, this commonly means defining as much as 10 regions, thus requiring the creation of as much as 10 search trees per simulation, which is both processing time- and memory intensive.
Finally, the SNESIM method as published attempts to generate the search tree or trees (in the case of multiple simulation regions), and run a simulation of the property field in a single computational workflow. This means in practice that if the input parameters that control the construction of the search tree are changed, a separate simulation must be performed for each change in the search tree. This precludes the use of SNESIM in an automated workflow, for example, to assess sensitivity to variations in training image input parameters, or a workflow where multiple training images are for a given set of conditioning data or varying amounts of conditioning data.